
clear
set more off
cap log close

***************************************************************************************************
* Program: figureS2.do
* Purpose: Produce Figure S2 (Smartphone Counts vs. Employment Level (QCEW))
***************************************************************************************************

************************************************************************
** set macros (directory)
************************************************************************
global prison_dir "~/Dropbox/Prison_Covid/PNAS_Nexus_Replication"
global raw_data "${prison_dir}/raw_data"
global intermediate_data "${prison_dir}/intermediate_data"
global code "${prison_dir}/code"
global analysis_data "${prison_dir}/analysis_data"
global output "${prison_dir}/output"

************************************************************************
** Main Code
************************************************************************

use "${intermediate_data}/prison_phones_count_June.dta", clear

* summarize at the county level 
gegen phones_count_cty = total(phones_count), by(COUNTYFIPS)

* collapse at county level 
sort COUNTYFIPS 
drop if COUNTYFIPS == COUNTYFIPS[_n-1]
keep COUNTYFIPS phones_count_cty 

* merge w/ BLS data
rename COUNTYFIPS area_fips
merge 1:1 area_fips using "${intermediate_data}/CA_92214_QCEW_BLS_062020.dta"
drop if _merge == 2

* calculate percentage
gen percentage = phones_count_cty / emplvl_062020

* Tobit / Heckman Fit
tobit emplvl_062020 phones_count_cty, ll(979) noconstant
predict fitted_tobit

replace emplvl_062020 = . if emplvl_062020 == 0
heckman emplvl_062020 phones_count_cty, select(qtrly_estabs) noconstant
predict fitted_heckman

* Scatter plot
tw (scatter emplvl_062020  phones_count_cty , xtitle("Smartphones Count") ytitle("Employment Level (QCEW)") ///
	text(7300 1055 "Heckman: y = 3.266x") text(7700 1020 "Tobit: y = 4.830x") ) ///
	(line fitted_tobit phones_count_cty, lpattern(dash) ) ///
	(line fitted_heckman phones_count_cty, lpattern(dot) legend(order(2 "Tobit fit" 3 "Heckman fit")))  
graph export "${output}/figureS2.png", replace 	
